### Script for examining why isolated PVA populations decline to extinction ###

# Load functions from PVA.R ##

#mother modified to show sex ratio
mother.mod<-function(n=2000, spX=10, spY=10, alpha=-1, fsurv1=0.4, fsurv2=0.1, msurv=0.1, init.b=0, init.b.var=3, h=0.3, gens=40, beta=100, fec=4, prob.d=0.1, sel.time=100, plot=FALSE){
	pop<-init.inds(n, spX, spY, init.b, init.b.var, h)
	n.list<-neighbours.init(spX, spY)
	sex.ratio<-c()
	sel<-FALSE
	for (g in 2:gens){
		pop<-repro(pop, spX, spY, beta, fec, init.b.var, h)
		if (g>sel.time) sel<-TRUE
		pop<-age(pop, selection=sel, alpha=alpha, fsurv1=fsurv1, fsurv2=fsurv2, msurv=msurv,spX=spX, spY=spY, beta=beta)
		pop<-disperse(pop, n.list, prob.d, spX)
		if (length(pop[,1])==0) {
			sex.ratio<-c(sex.ratio, NA)
			print(paste("The population went extinct at generation ", g, ":-("))
			break	
		}
		sex.ratio<-c(sex.ratio, sum(pop[,"S"])/length(pop[,1]))
		if (plot==T) plotter(pop, popsize, spX, spY, sel.time)
	}
	list(pop, sex.ratio)	
}

###################################################################################
## Model Runs ##
# fitted values
load("/Users/ben/Documents/Papers/EvoPVA/Quoll model R/Model inference/Posteriors Summary.RData")
ps<-out
rm(out)


	#Set parameters
	alpha.s<-ps["tf.alpha","Median"]
	fsurv1.s<-ps["tf.fsurv1","Median"]
	fsurv2.s<-ps["tf.fsurv2","Median"]
	msurv.s<-ps["tf.msurv","Median"]
	beta.s<-ps["tf.beta","Median"]
	fec.s<-ps["tf.fec","Median"]
	
	temp<-mother.mod(n=30, spX=1, spY=1, gens=100, alpha=alpha.s, fsurv1=fsurv1.s, fsurv2=fsurv2.s, msurv=msurv.s, beta=beta.s*1000, fec=fec.s, prob.d=0)
